Eddie module stream discharge from handout

This script examines the discharge in Neversink Cr.

This script uses the USGS data retrieval package to download data fro Neversink Cr. This is more up to date than the excel sheet.

Install libraries.

These are the libraries that are used in the script. These are installed only once and then can be commented out so they do not run in future execuations of the script.

# https://github.com/USGS-R/dataRetrieval

# general R packages to install if you dont have them
# remove the "#" if you need to run the installations

# install.packages("devtools") # essential in installing other thigs
# install.packages("tidyverse") # installs a lot of things and ggplot
# install.packages("scales") # allows great scale formatting on ggplot
# install.packages("lubridate") # makes working with dates easier
# install.packages("janitor") # clean names of columns and other things
# install.packages("readxl") # read in excel files
# install.packages("plotly") # interactive plots
# install.packages("skimr")
# install.packages("patchwork") # multiple plots

# specific to this module
# install.packages("dataRetrieval") # USGS Data Retreiveal Method

Load Libraries

This code loads the libraries for working with the data

# load these every time  
library(tidyverse) # the tidyverse library with many sub libraries
library(scales) # working with ggplot scales 
library(lubridate) # working with dates
library(janitor) # cleans up excel imports
library(readxl) # read in excel files
library(skimr) # does summary stats and other things
library(plotly) # makes interactive plots
library(patchwork) # combines many plots together

# # load the USGS package
library(dataRetrieval) # utilities for downloading and working with USGS data

Source of Neversink data

This is the data that is in the excel sheet that is used in this script.
https://waterdata.usgs.gov/ny/nwis/uv/?site_no=01435000&PARAmeter_cd=00065,00060,63160

This is all of the data that is available for the site.
https://waterdata.usgs.gov/nwis/inventory/?site_no=01435000&agency_cd=USGS

Read in USGS Data

This allows you to read in the daily mean data for the USGS site.
Note Enter the following information and you can change this to get different sites or time ranges:
site_no: this is the site that you want to import
par_code: these are the paramter code numbers - 00060 is discharge
You can retreive other variables using:
c(“00060”, “00065”)
00010 water temp
00060 discharge
00065 gage height
00095 specific conductivity if available

start.date: date to start on
end.date: date to end on

# The site is Neversink and the site number is 01435000
# USGS 01435000 NEVERSINK RIVER NEAR CLARYVILLE NY

# enter these variables
site_no <-          "01435000"
par_code <-         "00060"
start.date <-      "1800-01-01"
end.date <-        "2050-01-01"

neversink_daily.df <- readNWISdv(siteNumbers = site_no,
                        parameterCd = par_code,
                        startDate = start.date,
                        endDate = end.date)
  
# rename the columns
neversink_daily.df <- renameNWISColumns(neversink_daily.df)

Rename columns

This just renames the columns to lower case and variable names with units

neversink_daily.df <- neversink_daily.df %>%
  select(agency_cd, 
         site_no, 
         date = Date, 
         discharge_cfs = Flow) 

Create a month variable

The module would like to look at flows during differnet months. This creates a month column that can then be used to filter out and plot only the respective months.

neversink_daily.df <- neversink_daily.df %>%
  mutate(month = month(date))

Plot all discharge data

This is a plot of all of the data that is on the USGS website for daily mean values.

# saves the plot as q.plot so can be used later on   
q.plot <- neversink_daily.df %>% # use this dataframe
  ggplot(aes(date, discharge_cfs)) + # plot discharge_cfs versus date
  geom_point(size=0.1) + # add points that are small
  geom_line() # add a line on top of points
q.plot # call the plot to display

Plot February data

This creates a plot of just February data. It can do this using pipes rather than creating subdataframes which is a key advantage of tidyverse.

# saves as the feb.plot using the neversink_daily.df data
feb.plot <- neversink_daily.df %>% 
  filter(month == 2) %>% # this filters out only February
  ggplot(aes(date, discharge_cfs)) + # this plots discharge versus date
  geom_point(size = 0.1) + # adds small points
  geom_line() + # adds a line
  labs(y="Mean Monthly Discharge m^3/sec", x="Date") # adds clean axes labels 
feb.plot

NA

Interactive plot of feb data

This is an interactive plot of the february data. You can zoom in and out and select particular pints to see the values.

ggplotly(feb.plot)

August data

Plot of August data only using eh above formate. The only thing to change is the value for the month variable. This could be improved by showing the monthly mean of data in the same format below adding a a stat_summary command to show monthly mean as well.

aug.plot <- neversink_daily.df %>%
  filter(month == 8) %>%
  ggplot(aes(date, discharge_cfs)) +
  geom_point() +
  geom_line() +
  labs(y="Mean Monthly Discharge m^3/sec", x="Date")
aug.plot

Compare the February and August plots

We can use patchwork to put both plots on the same page to compare them

feb_aug.plot <- feb.plot +
                aug.plot +
                plot_layout(ncol = 2)
feb_aug.plot

Add trend lines of discharge over time

Because the data is over such a wide range I have zoomed in. There are two ways to zoom - the scale_x_continuous(limits = c(0,300)) will delete data that is not in this range. The coord_cartesian(ylim=c(0,300)) will zoom in without deleting data.

feb_aug.plot <- feb.plot + 
                geom_smooth(method="lm") + 
                        coord_cartesian(ylim = c(0,300)) +
                aug.plot + geom_smooth(method="lm") + 
                        coord_cartesian(ylim = c(0,300)) +
                plot_layout(ncol = 2)
feb_aug.plot

Regression statistics of February data

feb.model <- neversink_daily.df %>% 
  filter(month == 2) %>%
  lm(discharge_cfs ~ date, data=.)
summary(feb.model)

Call:
lm(formula = discharge_cfs ~ date, data = .)

Residuals:
   Min     1Q Median     3Q    Max 
-129.4  -86.5  -55.4    0.3 5925.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.585e+02  5.282e+00  30.007   <2e-16 ***
date        1.393e-03  5.713e-04   2.439   0.0148 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 233.5 on 2258 degrees of freedom
Multiple R-squared:  0.002627,  Adjusted R-squared:  0.002186 
F-statistic: 5.948 on 1 and 2258 DF,  p-value: 0.01481

Regression statistics of August data

aug.model <- neversink_daily.df %>% 
  filter(month == 8) %>%
  lm(discharge_cfs ~ date, data=.)
summary(aug.model)

Call:
lm(formula = discharge_cfs ~ date, data = .)

Residuals:
   Min     1Q Median     3Q    Max 
 -91.5  -61.5  -36.9   -2.6 7104.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 8.495e+01  4.794e+00  17.720  < 2e-16 ***
date        1.987e-03  5.152e-04   3.856 0.000118 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 220.2 on 2477 degrees of freedom
Multiple R-squared:  0.005966,  Adjusted R-squared:  0.005564 
F-statistic: 14.87 on 1 and 2477 DF,  p-value: 0.0001184